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Abstract 

We numerically investigate gauge field instabilities in anisotropic SU(2) plasmas using 
weak field initial conditions. The growth of unstable modes is stopped by non-abelian 
effects for moderate anisotropy. If we increase the anisotropy the growth continues be- 
yond the non-abelian saturation bound. We find strong indications that the continued 
growth is not due to over-saturation of infrared field modes, but instead due to very 
rapid growth of high momentum modes which are not unstable in the weak field limit. 
The saturation amplitude strongly depends on the initial conditions. For strong initial 
fields we do not observe the sustained growth. 
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1 Introduction 



It is an interesting and still open question to what degree the medium created in 
high energy heavy ion collisions reaches local thermal equilibrium before it falls apart. 
Results from the RHIC experiments are hinting towards fast thermalization [1]. For 
sufficiently large collision energy the relevant running coupling is small and this problem 
can be addressed theoretically in a controlled way using perturbative QCD. Even if a 
full analytic calculation is not possible one should at least be able to obtain parametric 
estimates for the thermalization time and the achieved temperature. Remarkably, the 
solution to this problem has not been found yet. 

Due to the non-isotropic expansion the momentum distribution of the produced par- 
tons becomes anisotropic^. If the expansion is mostly 1-dimensional along the collision 
axis, the typical longitudinal momenta become much smaller than the transverse mo- 
menta. Anisotropic momentum distributions cause so called plasma^ instabilities, i.e., 
certain long wave gauge field modes grow exponentially so long as their amplitudes are 
sufficiently small. This is a collective phenomenon which is not visible in the kinetic 
equation approach used in [3, 4, 5, 6]. It has been argued that this effect, which is well 
known in plasma physics, will speed up equilibration in heavy ion collisions since the 
unstable modes tend to make the momentum distributions more isotropic [7]. 

There are important qualitative difi^erences between QED and QCD plasma instabil- 
ities [10]. In both cases the growth of unstable modes is stopped by non-linear effects. 
In QED this happens when the amplitude of the unstable modes has become so large 
that they deflect a particle momentum by a large angle within a distance of one wave- 
length. This corresponds to gauge fleld amplitudes A of order p/e where p is a typical 
particle momentum, henceforth called " hard" . When the flelds become this large they 
have a dramatic effect on the plasma particles since they instantaneously make the 
momentum distribution isotropic. In QCD the gauge fields are self-interacting, and 
the linear approximation already breaks down at much smaller amplitudes A k/g 
where A; <^ p is a characteristic wave vector of an unstable gauge field mode. A crucial 
question is whether these non-linearities stop the growth of instabilities. In Ref. [8] it 
was suggested that gluon self-interactions may not saturate the instabilities because 
the system can "abelianize" so that the unstable modes can grow until they hit the 

^For a nice illustration see Fig. 1 of Ref. [16]. 

"^Hcrc "plasma" refers to a system of quarks and gluons which is not necessarily in thermal equi- 
librium, while sometimes the term "quark-gluon-plasma" is reserved for thermalized or almost ther- 
malized systems. 
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abelian saturation bound A<p/g. The distribution of hard gluons would then quickly 
become isotropic, and it has been argued [9] that this is sufficient for a hydrodynamic 
description to be applicable even if there is no local thermal equilibrium. 

The question how plasma instabilities in QCD saturate is thus an important one. It 

can be addressed most cleanly by neglecting both the expansion of the system and the 
back reaction on the particle momenta. This is sensible because in the weak coupling 
limit^ the expansion is slow compared to the dynamics of the unstable modes and 
because there is a large scale separation of particle momenta p and the wave- vectors of 
unstable modes k. 

Because the amplitudes of the unstable field modes become large, we are dealing 
with a non-linear problem and we cannot compute their time evolution perturbativcly. 
So far our qualitative understanding is very limited and one has to rely on lattice 
simulations. These are possible due to the large occupation numbers which allows one 
to use the classical field approximation for the infrared fields. In lattice simulations with 
fields depending only on t and z it was indeed observed [15] that the fields continue to 
grow rapidly in the non-linear regime. However, 3-1-1 dimensional simulations [11, 12] 
indicate that the instabilities are saturated by non-abelian interactions which would 
mean that their effect is less dramatic than suggested in Ref. [9]^. In [14] it was shown 
that even then the thermalization process is affected by plasma instabilities, because 
the broadening of longitudinal momenta of the particles caused by the unstable modes 
is more efficient than due to elastic scattering [6] . 

Most lattice simulations have so far been restricted to moderate anisotropics. In 
the present article we report on the evolution of instabilities in strongly anisotropic 
systems. In Sec. 2 we describe the equations and the approximations we use to solve 
them. The results are discussed in Sec. 4. In Ref. [20] strongly anisotropic plasmas 
have been considered in a kinematics and with approximations which are quite different 
from ours. 

^More precisely, one has to consider not only weak gauge coupling but also sufficiently large times 
where the system is sufficiently dilute so that the very notion of particles is applicable. In this regime 
the expansion rate is parametrically small compared to the time scale relevant to the instabilities 
[6, 10]. 

^For a recent discussion of the role of dimensionality see [13]. 
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2 The setup 



Our starting point is the non-abelian Vlasov equations [18, 19] 




(1) 



(2) 



These are classical equations of motion for SU(2) gauge fields interacting with 

particle degrees of freedom. The average distribution of the particles f{p) > is a 
gauge singlet, and the leading charged particle density fluctuations are described by 
adjoint representation distribution functions f"'{x,p). The particles are moving with 
the speed of light, thus, the 3-velocity is v = and (v^) is defined as {l,v). 

We neglect the back reaction of the soft gauge field on / and also the expansion, 
so we take f{p) to be space and time independent. Neglecting the x^-dependence of / 
is justified as long as the expansion rate of the system is small compared to the growth 
rate of the unstable modes we are interested in. In an isotropic plasma / only depends 
on \p\; here we consider the anisotropic case, but we assume that / is invariant when 
p is reflected or rotated around the z-axis. 

Our equations describe high momentum modes which are treated as classical colored 
particles and soft gluons which are treated as classical fields. In order for the classical 
particle approximation to be valid the wave vectors of the fields have to be much smaller 
than the momenta of the particles. The classical field approximation is valid because 
we will be dealing with large occupation number (large amplitude) gluon fields. The 
expansion of the system has been neglected because at weak coupling the expansion 
rate is much smaller than the rate at which the soft gluons evolve. Furthermore, the 
back-reaction of the soft fields on the momentum distribution has been neglected here 
('hard loop approximation'). 

The |p|-dependence of /" is irrelevant for determining the gluon field dynamics. One 
only needs the integral 



oo 




(3) 







Integrating (2) over |p| we obtain 




(4) 
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(v ■ Dwy = v^'F^y (5) 

with 



oo 





For isotropic / one would have u — rn^v, and (2) would the usual hard thermal loop 
equation of motion. For an anisotropic plasma u will not simply be proportional to v. 
Since we assume / to be parity even, u is parity odd. 
As in [21] we expand W{x, v) in spherical harmonics, 

Wix, v)^J2Yl Wim{x)YiM (7) 

1=0 m=-l 

with a finite Z-cutoff Lmax- This turns Eqs. (4), (5) into classical equations for fields 
living in 3+1 dimensions. Similarly we expand / in spherical harmonics and we assume 
that it only depends on and p^. Then 

Kp) - E M\p\)YiAv) (8) 

1=0 

where the sum runs over even I only. In general the Z-cutoff Lasym would be infinite, but 
in practice we must choose parametrizations with finite Lasym since the equations of mo- 
tion limit Lasym < L^ax- When we increase Lgsym it becomes possible to describe more 
anisotropic distributions, but at the same time L^ax and correspondingly memory- 
and cpu-time requirements of the simulations are increased (roughly proportionally to 

7-2 \ 
max/ ■ 

The equations of motion in terms of Wim in temporal gauge Aq — become 

doWim + q^^,^,D^Wi,„, = Foy^ + 2F,^ (9) 
doF^' + DkF'^'' = vl^W^m. (10) 



Gauss law reads 



Here = — is the canonical momentum of the gauge field A''. 



DiF'^ = -^Woo. (11) 

V47r 
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The coefficients Q^^/^/ may be found in Appendix A of Ref. [21]. Tfie otfier coeffi- 
cients are 

<,^l<Myzy, (12) 

We now define 

oo 

mf^AV^g^J 0^J,{p) (13) 



For an isotropic system ml equals the Debye mass squared. We want / to be positive 
which gives the condition '^i'mfYiQ{v) > (Albeit we shall violate this condition 
slightly.). 

The only non-vanishing w-coefiicients in Eq. (9) are 



^ /(/ + 1) ( ml^ ^Hi \ 



uE = V^(7Tl)m,2 (17) 



u^i^_, = < = = (18) 

We study the behavior of the system using both weakly and strongly anisotropic 
distributions. A measure of the anisotropy is 

ry^ ^ Hvl)/{v') , (19) 

which equals 1 for symmetric and for completely planar distribution. 

For each Lasym the distribution is parameterized by the coefficients mf, with / = 
0,2,..., Lasym- The values of mf are chosen so that the anisotropy of the resulting 
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Figure 1: Anisotropic hard particle distributions used in this work, together with the 
distribution used by Arnold, Moore and Yaffe [11]. The distributions are plotted so 
that the relative number of particles moving to direction v is proportional to the length 
of the radial vector from the center of the plot. For each Lasym we tried to maximally 
localize the distribution in the xy-plane. The distributions are normalized to equal 
area for readability. 

distribution is approximately maximized. The reason for this choice is that for a given 
anisotropy, we take Lasym as small as possible, also minimizing the required Lmax and 
hence computational requirements. 

For Lasym = 2 and 4 the tuning of the parameters is easy enough to do by hand, but 
for Lasym = 14 and 28 we use a 1-parameter fitting procedure: / (9) is fitted to a narrow 
Gaussian function centered at = 7r/2. The width of the Gaussian is adjusted to be 
as small as possible while still giving a good fit; if the width of the Gaussian is too 
small the fitted function will have large oscillations over whole 6'-range. The quality of 
the fit is justified by eye. This procedure is sufficient for our purposes: the goal is to 
find one good enough parametrization for the asymmetry, and no attempt is made to 
maximize the asymmetry for any given Lasym- The resulting parameters are given in 
table 1. 

This process gives distributions where the power is strongly concentrated around 
6 = 71/2 ± A6, where is the maximum resolution power of the Y/o-^xpansion when 
i < -^asym, that is A^^ ~ 7r/Lasym- Thus, wheu plotted on cartesian coordinates, the 
distribution has well-defined "lobes" centered around direction 9 = n/2, i.e. along the 
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Table 1: The parameters used in simulations. L^syxa — 2 and 4 correspond to weak 
asymmetry, Lasym = 14 and 28 to strong asymmetry. 

xy- plane, as shown in Fig. 1. For directions near ^ f« or tt, the distributions can 
become slightly negative; however, the magnitude of this effect is negligible. 

For small amplitudes the non-linear terms in the equations of motion can be ne- 
glected. Modes with different wave vectors do not mix, and the unstable modes grow 
exponentially at a rate which can be calculated analytically. The growth rate is shown 
in Fig. 2 as a function of the length of the wave vector of the unstable mode for different 
asymmetries. For each asymmetry A;* denotes the value of |fc| for which the growth 
rate is maximal. For the smallest to the largest anisotropy, the maximum growth rate 
increases by a factor of 5 and the width of the unstable mode distribution by a factor 
of 8. 

The linear equations of motion offer a straightforward method for investigating how 
large we need need to make in order to reproduce the continuum dynamics. In 
Fig. 3 we compare the growth rate at ~ oo with the rates at different finite values 
of ivmax for modes with k — kz. For weak anisotropy (Lasym ~ 2, left figure) one 
needs rather large values of Lmax ^ -^asym to reproduce the growth rate. The growth 
rate for strong anisotropy Lasym = 28 (right figure) can be reproduced already with 
-^max ^ -^asym- Indeed, for the asymmetries used in this study it appears that the finite 
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Figure 2: Soft field growth rate as a function of momentum k = kz for linearized 
equations of motion, for anisotropic hard mode distributions Lasym = 2, 4, 6, 14 and 
28 (see Fig. 1 and table 1). /c* is the wave number with the maximal growth rate. 
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Figure 3: Growth rate of magnetic energy for the linearized equations of motion with 
different Lmax cutoffs, shown for Lasym = 2 (left) and Lasym = 28 (right). 
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-^max effects are rouglily independent of the Lasym used, and we should obtain accurate 
results for Lmax ^16, of course provided that we keep Lmax > -^asym- In Sec. 5 we 
investigate the Lmax-dependence of the real simulations in detail. 

We note that while the rate can be solved analytically at finite Lmax, in Fig. 3 we 
actually measured the rate from numerical simulations using a linearized version of 
our simulation program. Thus, this measurement was also an important check of the 
correctness of the simulation program. 

3 Simulation program and parameters 

The equations of motion (9), (10) are discretized as described in Ref. [21]; we invite 
interested readers to check therein for the detailed implementation. We note that 
we implement the W^-fields in a "staggered" fashion: because the W equations of 
motion only have first order derivatives, a symmetric discretization decouples ly-fields 
at even and odd lattice sites from each other (i.e. at space-time sites where the integer 
coordinate rix + Uy + Uz + rit, where rit is the number of the evolution time step, is 
either even or odd.). Thus, we can delete the 1^- field at odd sites, saving memory and 
cpu-time.^ 

For the time update we use a time-symmetric staggered leapfrog as described in [21]. 
The only essential difference is the appearance of the last term in Eq. (9). In order 
to guarantee that the update remains invariant under time reversal we implement the 
update of the VF-fields in two stages, interleaving these with the gauge and electric 
field update steps. 

The time-step values we use are St = 0.05a and 0.1a, where a is the spatial lattice 
spacing. We shall discuss the lattice artifacts - finite a, finite volume, finite 6t, and 
finite L-^ax ~ in detail in Sec. 5; to summarize, all lattice effects appear to be well under 
control. 

We note that while all mf are dimensionful in the equations of motion, for fixed 
asymmetry the ratios mf/iriQ remain constant. Thus, every dimensionful quantity can 
be given in terms of the powers of single parameter, mg. In particular the lattice 
spacing is given as (amg). The gauge coupling constant g"^ can be completely absorbed 
in the equations of motion, making the results independent of the value of g^. 

^This procedure also deletes half of the unphysical doublers inherent in the M^-field spectrum. The 
reason these doublers appear is the same as for the notorious lattice fermion doublers, namely the 
first order derivatives. However, in our case the doublers are quite benign, as is discussed in [21]. 
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Table 2: The lattice spacings (in units of mo) and lattice sizes used in the weak initial 
field analysis for each value of the asymmetry. For several of the volumes there are 
more than one individual run. The Lmax-cutoff used is shown at the top of the columns. 
In addition, there are some some volumes with more than one Lmax-cutoff; these are 
indicated with a subscript (only for Lgsym = 14, 28). 

Our initial conditions are as follows: we initialize the electric field components E"'{x) 
to a small amplitude white noise, i.e. random Gaussian fluctuations, with vanishing 
initial A and Wim- We make an orthogonal projection of the E'-fields to a hypersurface 
satisfying Gauss' law, DiE'^ = (since Woo = 0). The evolution equations preserve 
Gauss' law. The electric fleld drives the gauge flcld A to a non-zero value very quickly, 
so that (-B^) ^ {E'^) before the exponential growth of the unstable modes becomes 
visible. The amplitude of the initial fluctuations is chosen small enough so that the 
equations of motion arc essentially linear during the initial stage. The growth of un- 
stable modes then drives the flelds to much larger values. 

The lattice spacings and sizes used in the analysis are shown in table 2. The simu- 
lations have been performed mostly using pc-clusters with infiniband interconnects. 
The simulations require unusually large amounts of memory (for lattice simulations); 
our largest simulations used 192 nodes, with a total memory requirement of around 
400GB. The simulations were performed at the Finnish IT Center for Science (CSC). 
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4 Results 



4.1 Energy densities 

As mentioned in Sec. 3, the initial condition we use is a white noise spectrum satisfying 
Gauss' law for the electric field, with A and W set to zero. In Fig. 4 we show the field 
evolution for weak anisotropy (Lmax = 2,4) starting from very small amplitude initial 
conditions for different values of the lattice spacing. We see quahtatively the same 
behavior as observed in Refs. [11, 12]. After some initial setthng down, the soft fields 
start to grow exponentially until they reach the non-abelian point A ^ hj g where non- 
linear terms in the equation of motion start playing a role. We find that this happens 
when the magnetic field squared approximately equals 

1 h'^ 

-B^ (20) 

2 non— abclian Acj'^' \ ' 

After that the growth slows down significantly and is no longer exponential. According 
to Arnold and Moore [16] this growth is due to cascading of energy from the originally 
unstable infrared modes to higher momentum ones. The amplitude of the initial fields 
was not specifically tuned to be equal for different lattice spacings; nevertheless, the 
gauge field evolution falls on a curve independent of the lattice spacing (as long as the 
volume is large enough, see Sec. 5). The origin of time t = has been adjusted in 
Fig. 4 in order for the growth phases to overlap. Thus, only differences of t have a 
physical meaning. 

For strong anisotropics we find a very different picture. In Fig. 5 we show our results 
for Lasym = 14 and 28. We clearly see the onset of non-linear effects at the magnetic 
field energy density around kl/{Ag^). There the growth ceases to be exponential and 
the dynamics becomes very complicated. The electric field grows very rapidly, and the 
electric field energy becomes as large as the magnetic one. Subsequently, however, the 
growth of energy continues at a large rate. It is not a purely exponential growth, but 
it is not much slower than the initial weak field growth. For Lasym = 28 the growth 
rate is roughly as large as in the weak field regime (mot < 40). 

At some value of the energy density the growth saturates. Furthermore, in contrast 
to the moderate asymmetry in Fig. 4, the electric and magnetic field energies reach an 
equal level at the end. In Fig. 5 we show the values where the growth finally saturates 
for different values of the lattice spacing a. We see that the saturation energy has a 
strong dependence on the lattice spacing, growing as amo is decreased. Therefore we 
can conclude that the saturation seen in Fig. 5 is caused by the lattice regularization. 
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Figure 4: Magnetic and electric field energy densities as a function of time for moderate 
anisotropy, measured from lattices with different lattice spacings a. The lattice sizes 
are the largest ones for each lattice spacing in table 2. 

In Fig. 6 we show the maximal magnetic energy density as a function of the lattice 
spacing. The maximal energy density appears to grow without bound with decreasing 
a with a power-like behavior. The magnetic field energy on the lattice is given by 
4/ {ag^) X]i<j(l ~ |Tr t/jj), where Uij is the ordered product of link variables around a 
spatial plaquette, 

Uij {x) = U, (x) Uj {x + al) U} (x + a] ) U]{x). (21) 

There is an absolute upper limit on the magnetic energy density, 24/(a^(7^) which is 
reached when Tr Uij = —2. This is a very particular fully ordered state; a more realistic 
limit is the completely random state where (Tr Uij) = and where the magnetic energy 
density reaches the limit 12/{a'^g'^). Energies above this limit are shown in Fig. 6 as a 
shaded region. 

We observe that our maximal field energies do not quite reach the maximum energy 
limit. Instead, the saturation energy density appears to diverge in the continuum limit 
with a different power of a. If we fit a power law behavior to the saturation energy 
density at both asymmetries, we obtain the results -^saturation oc (amo)~^'^ for L^ax = 14 
and (amo)"^'^ for L^ax = 28. Because we do not have proper statistical errors for the 
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Figure 5: Same as Fig. 4 but for stronger anisotropy. Now the growth of field energy 
appears to continue indefinitely and it is stopped only by lattice cutoff effects. For 
each lattice spacing we show the largest volume listed in table 2. 



data in Fig. 6, we cannot quote proper error bars for the fitted exponents. However, 
we can nevertheless make a rough estimate of them by performing jackknife analysis 
in terms of the individual simulation points, obtaining an error bar ±0.2 for both 
exponents. It is worth noting that the exponent in the L 



asym 



28 case is close to —3, 
the exponent given by the thermal distribution with a lattice cutoff. 

This analysis shows that there appears to be no saturation of the energy density if 
the lattice spacing is removed. This is very different behavior from the one that was 
observed in the 3+1 dimensional simulations of Refs. [11, 12]. 

Let us now discuss possible reasons for this behavior. When the anisotropy is mild, 
the unstable modes have momenta of order itiq. However, for strong anisotropy there 
are unstable modes with \k\_\ <mo but with longitudinal momentum \kz\ all the way 
up to /cmax, wherc 
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(22) 



and 1] is the measure of anisotropy introduced in eq. (19). In [23] it was argued that 
the magnetic field squared of these modes cannot become larger than ~ tuq/ {g^rf')- 
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Figure 6: Maximal magnetic field energy density as a function of the lattice spacing 
for Lasym = 14 and 28, and for all values of Lmax used. The shaded region is above the 
maximum magnetic field energy density, given by a completely random lattice gauge 
system. The dashed lines are power-law fits to the two asymmetries, with the results 
(amo)"2-^ (Lasym = 14) and (amo)"^-^ (La^ym = 28). 

The energy density at saturation in a strongly anisotropic plasma would then be en- 
hanced by a factor 1/rf compared to the case of moderate anisotropy. However, this 
enhancement factor is only about 16 for Lasym = 14 and about 67 for Lasym = 28, while 
we see the energy density in Fig. 5 growing by many orders of magnitude larger than in 
the case of weak anisotropy. Therefore it is not a (quasi-) exponential growth of modes 
who's equations of motion are almost linear which could explain the behavior seen in 
Fig. 5. Thus the continued growth must be an effect which is essentially non-linear. 

There appear to be (at least) two scenarios for the physics behind the continued 
growth. The first is that the unstable modes grow to occupancy much larger than 
l/g"^ as suggested in Ref. [8]. Another possibility is that the energy goes into the high 
momentum modes, rather than into the modes which are unstable in the weak field 
regime. We shall try to distinguish between these outcomes by measuring quantities 
which are sensitive to the momentum spectrum of the gauge fields: gauge fixing and 
direct Fourier transformation, gauge invariant operators and gauge invariant cooling. 
These all indicate that the energy indeed gets dumped to the UV, and there is no 
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growth of the IR modes much beyond the non-abehan point. 
4.2 Coulomb gauge occupation numbers 

For free gluon fields the concept of occupation numbers /s(fc)^ is unambiguous. It can 
be calculated from the gauge field by fixing to Coulomb gauge using the expression 

(23) 

where A^dof denotes the number of color/spin degrees of freedom. For refiection invariant 
field configurations the interference term of A and E vanishes. For free fields the two 
remaining terms give equal results when they are averaged over time. Thus, assuming 
refiection invariance, one can compute the occupancy either from A or from E only, 
and in this work we use the former case. The distributions shown here are averaged 
over all directions of fc, 

/(^) ^ / ^/s(fc) (24) 

If the gluon field amplitudes are large and/or the gluons are interacting with the 
particles, there is no occupation number in the strict sense. Nevertheless one expects 
that (23) still gives a good estimate of the power in one field mode. However, fixing 
the gauge for large fields in a non-abelian theory is dangerous due to Gribov copies of 
near vacuum configurations of the high momentum modes. We make three consistency 
checks of the gauge fixed spectrum by comparing with gauge invariant measurements: 
the total energy in the gauge fixed spectrum, measurement of the average (fc^), and 
comparing the spectrum with gauge invariant cooling. These will be discussed below. 

The occupation numbers as a function of time are shown in Fig. 7, for strong (Lgsym = 
28) and moderate anisotropy (I/asym = 4). The curves show the spectrum measured 
at constant evolution time intervals. Early times are at the bottom; the initial white 
noise £?-field implies a spectrum /(/c) ~ l/Zc. 

Let us first consider the case of strong anisotropy. At early times one sees a rapid 
growth of the infrared modes which is the fastest at k = k^. The dashed curve is at the 
time at which non-linear efi^ects become important. In Fig. 5 this time is marked with a 
vertical dotted line. As this is happening the active mode spectrum widens very rapidly. 

^Wc use the subscript s to distinguish the occupation number of the classical (soft) fields from the 
occupation number of hard gluons which are described by the VF-field. 
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Figure 7: Coulomb gauge power spectrum (occupation number) as a function of time 
for strong (Lasym = 28, left) and weak (Lasym = 4, right) anisotropy. The power spectra 
are plotted at equal intervals of At = 3.6/mo for Lasym = 28 and At = 16.4/mo for 



At later time times the amplitude of the k ~ k^^-modes does not grow any longer, but 
the ultraviolet end of the spectrum grows extremely rapidly - in fact the occupation 
number at higher k grows faster than the original growth rate at k^, as can be observed 
from the large gaps between the lines in Fig. 7. The final spectrum is shown with a thick 
line, and its shape fits f{k) ~ l/Zc quite well, consistent with a thermal distribution. 
However, a more detailed inspection of the spectrum shows that the growth of the 
energy stops before this is reached: the growth stops when the occupation numbers 
near the lattice cutoff k/nio = 7r/(moa) ~ 10.5 become appreciable (> 0.05). After this 
the distribution just settles towards the thermal one, without increase in energy. 

The situation at modest anisotropy (Fig. 7 right) looks quite similar at the beginning. 
However, in this case the growth in the UV part of the spectrum stops soon after the 
non-abelian point is reached. The mode spectrum remains dominated by the IR modes, 
and the total energy grows only approximately linearly with time. 

As a check that the occupation number reflects the true distribution of energy over 
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Figure 8: Energy density computed from the Coulomb gauge power spectrum, com- 
pared with the true energy density in the magnetic field, for Lasym = 28, niQa = 0.3. 

the different modes we compute the total field energy density e 

and compare it to the (gauge invariant) direct measurement of the energy from the 
lattice. The result is shown in Fig. 8. In the weak field regime our measured / slightly 
over-estimates the energy density. One has to keep in mind that even for very small 
amplitudes the gauge fields are not free, but are coupled to the PF-fields, so that the 
two curves need not coincide exactly. At large fields the discrepancy is shghtly bigger 
and / yields a slightly too large result. However, the overall disagreement is within a 
factor of 1.4. 

4.3 Average \k\ from gauge invariant operators 

The Coulomb gauge occupation numbers strongly indicate that the continued growth 
seen above is due to population of high momentum modes. However, one may be 
concerned about gauge artifacts, because strong fields could produce fake high mo- 
mentum occupancy. Therefore, in order to be certain about our conclusion regarding 
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Figure 9: Average as a function of time, measured from the gauge fixed occupation 
numbers /(/c), and from the gauge invariant operator, Eq. (26), for the Lasym = 28, 
mofl = 0.3 -simulation shown in Figs. 7 and 8. The shaded region is the time interval 
when the non-linear growth of energy is occurring. 

the high momentum occupation, it is mandatory to investigate this result also using 
gauge invariant measurements. A measure for the typical momentum squared of the 
color-magnetic fields is 



In QCD there is also the commutator [Aj, i?^] contributing to (fc^). So it appears 
that large (fe^) does not necessarily imply that the typical of the magnetic field 
is large. However, in the 1-dimensional simulations [15] where the unstable modes 
grow indefinitely, the commutator terms were found to remain small in accordance 
with the abelianization picture of Ref. [8]. Thus we expect our (A;^) to be a good 
measure of the momentum of the modes. Note in particular that the commutator term 




(26) 



In electrodynamics this would equal 




QED — 



J k^\B{k)\^d^k 
f\B{k)\^d^k 
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is parametrically of the same size as the gradient term when non-hnear effects start 
playing a role. 

In Fig. 9 we show (fc^) as a function of time, both computed from the gauge invariant 
object (26) and from the Coulomb gauge occupation numbers. At early times {t < 
12/ mo), when the fields arc very weak, (fe^) is large because it is dominated by UV 
modes due to our white noise initial conditions. As soon as the unstable modes start 
growing they give the dominant contribution to (fc^) which is then of order kl- The 
two curves do not coincide which is not surprising since even for free fields they would 
in general not be identical. Once one is in the non-linear regime, the average fc^ 
increases rapidly. This is a clear signal of a rapid transfer of energy to high momentum 
field modes. When the lattice cutoff starts having an influence on the time evolution 
(^^0 ^45), the two curves start to deviate strongly. 

4.4 Cooling 

Another gauge invariant method for obtaining information about the gauge field spec- 
trum at a given physical time is to take the gauge field configuration at that time and 
let it evolve in the (unphysical) cooling 'time' r using the equation of motion 

drAi = DjFji (27) 

This reduces the gauge field energy monotonously. For weak fields the Fourier compo- 
nents in Coulomb gauge evolve like 

Ai{T, k) = exp(-rfc2)A(0, k). (28) 

Thus, the cooling has the largest effect on the high momentum modes and they are 
depleted first. Results for the coohng time dependence of the magnetic field energy are 
shown in Fig. 10 (full lines), measured at intervals At — 3.6/mo during the evolution 

of a system with strong anisotropy (Lasym = 28). 

For free fields with a thermal spectrum Eq. (28) gives the result Energy ~ r~^/^ for 
large enough r. This behavior is clearly visible at early time cooling curves, the bottom 
curves in Fig. 10.^ When we are in the linear regime where the unstable modes grow 
exponentially, practically all of the the energy is in the infrared modes, and the cooling 
takes more time to have any effect on the total energy. This is visible as horizontal 

^Our initial condition was small amplitude white noise for E, which is thermal by itself. This 
rapidly populates A-modes to an approximately thermal distribution. 



19 



L =28 am =0.1 V = 240^ At = 4.0/m, 

asym ( 



1 1 — I I I I I 1 1 1 1 — I I I I 1 1 




cooling time T / a 

Figure 10: Magnetic field energy (solid lines) during the cooling of the field configu- 
rations. The dashed lines are obtained from the 'cooled' Coulomb gauge occupation 
numbers, Eq. (29). The different curves are for physical times in intervals of 4/mo, 
with time increasing from bottom to top. The final curves for both cases are shown 
with thicker lines. 

lines in the middle part of the cooling plot. When the cooling time reaches r ~ 

the energy starts to decrease rapidly and the cooling curves develop a smooth shoulder. 

The results from the gauge invariant cooling can be directly compared with the 
Coulomb gauge fixed field mode spectrum. Because f{\k\) oc | A(fc)p, we obtain 'cooled 
occupation numbers' from 

Uoi{k,T)^e-''''^f{k). (29) 

From this we can calculate the corresponding energy density as a function of r. These 
are plotted in Fig. 10 with dashed lines. We observe that these match the gauge 
invariant cooling curves perfectly at initial times where the field amplitudes are small. 

However, at around t = 28/mo (7th curve from the bottom, see also Fig. 5) the 
system enters the non-linear evolution domain and the two curves start to separate. 
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Figure 11: Time evolution of magnetic field energy for different choices of the initial 
field amplitude. 



This is due to two effects: firstly, the gauge fixed occupation number calculates energy 
slightly incorrectly for large fields, especially in the infrared end of the spectrum. 
Secondly, for large amplitude fields the non-linear equations of motion make the cooling 
significantly less efficient in reducing the energy. Thus, in the linear approximation used 
in Eq. (29) the energy decreases much faster than with the gauge invariant cooling. This 
is clearly visible in Fig. 10. Nevertheless, the main features are the same: the 'shoulder' 
in the cooling curves moves towards smaller r, which implies that the ultraviolet modes 
become occupied. 

4.5 Non-weak field initial conditions 

So far we have only considered very weak initial fields. With such initial conditions 
only modes which have k very close to the 2;-axis get substantially excited because this 
is where the growth rate is the largest. By the time the equations of motion become 
non-linear, the field's momentum distribution is almost 1-dimensional. It does not 
mean, however, that our results are just what has been observed in 1 + 1 dimensional 
simulations [15], where the growth continues beyond the non-abelian saturation limit. 
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This is because in our 3 + 1 dimensional simulations with moderate anisotropy the 
growth saturates even for very weak field initial conditions (cf. Fig. 4 and Ref. [11]). 

Let us now consider larger initial fields (Fig. 11). In this case we use only the strong 
anisotropy lattices, Lasym = 28, and moa — 0.3. The electric fields are now initialized 
with an infrared-dominated spherically symmetric spectrum, {E{k)) oc cxp[— fc^/(0.6mo)^]. 
The initial electric field energy densities vary from 0.0032 /{g~'^mQ) to 14.1/((yf~^mo); 
from the electric field the energy is rapidly pumped in the magnetic fields, as is evident 
from the figure. Note that the initial momentum spectrum is dominated by modes 

Wc see that there is a very strong dependence on the size of the initial fields. If the 
fields start out near the non-abelian point (20) there is practically no growth. ^'^ This 
behavior is very different from the one observed in Ref. [17] where there is growth for 
large initial fields. We leave more detailed analysis for further study. 

5 Lattice artifacts 

When a new phenomenon is studied with lattice simulations, it is very important to 
quantify possible harmful discretization and finite volume effects. The very large range 
of scales makes this check especially crucial in this case. As we shall detail below, all 
lattice effects appear to be well under control. 

Lattice spacing: The effects caused by different lattice spacings a were already 
discussed above. As can be seen from Fig. 4, at weakly anisotropic hard mode dis- 
tributions the finite a effects are small - the small dispersion of the results is of the 
same magnitude than statistical deviations at fixed a. We made no effort to enforce 
physically equivalent initial conditions for different values of a. On the other hand, the 
finite lattice spacing effects were seen to be quite large for strong anisotropy, Fig. 5, 
due to the population of the ultraviolet modes. Even in this case there appears to be 
an universal lattice spacing independent evolution, which finite a simulations follow 
before they finally saturate. 

^°It should be noted that in this case the system is not dominated by single mode k « k^z; thus, 
the 'non-abelian limit' for energy density does not describe the properties of the system as well as 
before. Nevertheless, we keep this quantity for comparison. 
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Figure 12: The growth in magnetic energy for Lasym = 14, mQa = 0.55 runs using 
different volumes. The 3 largest volume curves are practically on top of each other. 

Finite volume: If the volume is too small, it can effectively lower the dimensionality 
of the system. Indeed, too small volume can cause too much growth. In Fig. 12 we 
show the evolution using 4 different volumes for Lasym = 14, mga = 0.55 -case. Except 
for the smallest volume the curves fall on top of each other. (The statistical dispersion 
between the large volume runs is very small due to the smallness of the random initial 
fluctuations.) Thus, neither the exponential growth nor the final saturation can be 
due to the finite size of the system. In general, we require system sizes L > 5(27r/fc*), 
except for the very smallest lattice spacing. 

Finite L^nax- We have also studied the Lmax-dependence of the field growth. In 
Fig. 13 we show the magnetic energy density evolution for Lasym = 28, ruQa = 0.3 on 
a 128^ lattice, using L^ax = 32 and Lmax = 48. In this case we used identical initial 
conditions. As can be seen, the evolution is almost identical. (See also Fig. 3.) 

Finite 6t: In addition to finite lattice spacing a, in simulations of equations of 
motion one has to check the finite update time-step effects. In this work we used 
St = 0.1a, and checked the stability of the results against 6t = 0.05 simulations with 
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Figure 13: The evolution of the magnetic field energy density for Lasym = 28, mga = 0.3 
-lattices, using Lmax-cutoffs 32 and 48. The initial conditions were identical for the two 
runs. 

otherwise identical setup. The results are in practice indistinguishable, showing that 
our original 5t = 0.1a is sufficiently small. 

6 Summary and discussion 

We have studied the dynamics of infrared gauge fields in anisotropic SU(2) plasmas in 
the so called hard loop approximation, i.e., neglecting the backreaction of the infrared 
gauge fields on the phase space distribution of the high momentum partons. Starting 
from weak field initial conditions we find a behavior which appears to be qualitatively 
different from what was observed previously for weakly anisotropic plasmas. The field 
energy grows until non-linear effects start playing a role, which slow down the growth. 
But then the growth resumes and appears to continue without limit and it is only 
stopped by the lattice cutoff. For very strong anisotropy it is almost as fast as the 
initial exponential growth. This continued growth is different in nature from the linear 
growth found in weakly anisotropic plasmas. 

We have studied gauge fixed occupation, gauge invariant operators and cooling. All 
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methods indicate that there is a rapid transfer of energy to field modes which have 
\k\ ^ kuiax- These are modes which have no instabihties in the weak field regime. 

For the largest anisotropy we find a growth rate in the strong field regime which is 
approximately the same as in the weak field regime. The growth in total energy persists 
even though the magnitude of the soft gauge modes with |fc| ~ appears to remain 
constant. The mechanism of the energy transfer from the hard modes (VF-fields) to 
gauge field modes with |fc| ^ A;* remains unknown. 

We would like to point out that the earlier 3-dimensional simulations show an inter- 
esting structure which has not been discussed so far. After a weak field regime with 
exponential growth the system enters a phase where the fields become strong and non- 
linear effects become important. But then, after a brief pause, the fields again start to 
grow rapidly, almost as fast as during the initial exponential growth. Only after that 
there is finally a saturation and the subsequent linear growth. To reiterate, even in the 
weakly anisotropic case there appears to be a 2-stage structure in the saturation. It 
is conceivable that the behavior we observed is qualitatively similar. However, we find 
that this continued growth lasts much longer when we increase the anisotropy of the 
system, maybe forever. 
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